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ABSTRACT. We present some relations that allow the efficient approximate inversion of 
linear differential operators with rational function coefficients. We employ expansions 
in terms of a large class of orthogonal polynomial families, including all the classical 
orthogonal polynomials. These families obey a simple 3-term recurrence relation for 
differentiation, which implies that on an appropriately restricted domain the differenti- 
ation operator has a unique banded inverse. The inverse is an integration operator for 
the family, and it is simply the tridiagonal coefficient matrix for the recurrence. Since 
in these families convolution operators (i.e. matrix representations of multiplication by 
a function) are banded for polynomials, we are able to obtain a banded representation 
for linear differential operators with rational coefficients. This leads to a method of 
solution of initial or boundary value problems that, besides having an operation count 
that scales linearly with the order of truncation N } is computationally well conditioned. 
Among the applications considered is the use of rational maps for the resolution of sharp 
interior layers. 
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1. Introduction 

The solution of constant coefficient ordinary differential equations with periodic bound- 
ary conditions is especially simple in the Fourier spectral representation, since differen- 
tiation of a smooth function is replaced by multiplication of its Fourier coefficient vector 
by a diagonal matrix. An analogous property is shared by Hermite polynomial expan- 
sions in unbounded domains. Other spectral representations give, in general, almost full 
triangular differentiation matrices. However, for polynomial families such as the Cheby- 
shev and Legendre, the matrices representing some commonly occurring operators, such 
as the Laplace operator in various separable geometries, are known to be reducible to 
simple, banded form through the use of appropriate banded preconditioners ([10, Ch. 
10], [8], [16]). The origin of most of such simplifications is found in the fact that the 
matrix operator for integration in any of the classical orthogonal polynomial families is 
tridiagonal [7]. 

In this article we show how to exploit the properties of the operator of integration 
to arrive at efficient spectral algorithms for the approximate solution of a large class of 
ordinary differential equations of the form 

(1) Lu = ^2(m k (x)D k )u = f(x) , = 6), 

k—0 

subject to the constraints 

Tu = c. 

where m* are rational functions of x, D k denotes k - th order differentiation with respect 
to x, T is a linear functional of rank n, and c £ Rn. (Typically, the constraints are 
boundary or initial conditions, but this is not necessary.) 

The problem of approximating solutions of Ordinary or Partial Differential Equations 
(O. or P.D.E.) by spectral methods, known as Galerkin approximation, involves the 
projection onto the span of some appropriate set of basis functions, typically arising 
as the eigenfunctions of a singular Sturm- Liouville (SL) problem. The members of the 
basis may satisfy automatically the auxiliary conditions imposed on the problem, such 
as initial, boundary or more general conditions. Alternatively, these conditions may be 
imposed as constraints on the expansion coefficients, as in the Lanczos r-method [13]. 

It is well known [4] that the eigenfunctions of certain singular Sturm-Liouville problems 
allow the approximation of functions in C°° [a, 6] whose truncation error approaches zero 
faster than any negative power of the number of basis functions (modes) used in the 
approximation, as that number (order of truncation N) tends to oo . This phenomenon 
is usually referred to as ‘spectral accuracy 1 [10]. The accuracy of derivatives obtained by 
direct, term-by-term differentiation of such truncated expansions naturally deteriorates 
[4] but for low order derivatives and high enough order of truncation this deterioration 
is negligible, compared to the restrictions in accuracy introduced by typical difference 
approximations. Since results on the accuracy of spectral methods are well documented 
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in the literature, we shall limit ourselves to the discussion of certain formal properties of 
orthogonal polynomial families which allow algorithmic simplifications in their use. Facts 
about orthogonal polynomials that we shall need can be found in any of the standard 
references (e.g. [14], [17]). 

Throughout, we assume that we are working with a family of polynomials {Qfc}o° which 
are orthogonal and complete over the interval (a, b) (here a and/or b can be infinite) 
with respect to the nonnegative weight w(x). In the cases of interest, these are the 
eigenfunctions of a Sturm-Liouville problem 

(2) (p(x)Q' k )' + A k w{x)Q k = 0. 

Then the Q' k form an orthogonal family as well, with nonnegative weight p(x) which 
satisfies p(x ) — > 0 as x — *• o, b. In this paper we focus exclusively on the classical orthog- 
onal polynomials, i.e. the Jacobi, (special cases of which are the Chebyshev, Legendre 
and Gegenbauer polynomials) Laguerre and Hermite polynomials, which are the only 
polynomial solutions of Sturm-Liouville problems of the form (2) [12]. We will assume 
that the functions under consideration possess sufficient differentiability properties over 
(a, b) and can be expressed as a series involving the Q k - See [4] for a discussion of the 
convergence properties in the relevant function spaces. 

We introduce the spaces Q" by 

Q n m = span{Qk\m < k < n). 

Our method constructs an approximate particular solution of (1) in a subspace of codi- 
mension 7 T, (e.g. ( Qq *)■*■) such that when n-th order differentiation is restricted to this 
subspace it has a simple inverse. We also require that L be invertible when restricted to 
this subspace and that T has full rank when restricted to the space of solutions to the 
homogeneous problem, ((1) with / = 0). 

Of key importance for our purposes is the requirement that differentiation or its in- 
verse (‘integration’ in an appropriately restricted domain) must have banded form. For 
example, the first derivative operator in the Chebyshev representation, D has elements 

'0, i > j, 

1 _ I i < j, i + j even 

~ < j, 0 < i < j, i + j odd 

k i , i = 0, j odd 

Its inverse, when respective domains and ranges are appropriately restricted, is given 
by: 

0 0 0 0 ••• 0 \ 

2 0-1 0 

B= \ 0 1/2 0 - 1/2 0 . 

2 o o ••• ••• ••• 0 

0 0 0 1/Jfc 0 -1 Ik ••• / 
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Now, DB = /go o while BD = Jg°°. Clearly, D k B k = Iq«> as well. However, B k D k ^ I. 
If we apply &-fold differentiation to an arbitrary function, all information about the first 
k coefficients in its Chebyshev expansion is lost. If however we restrict the action of D k 
to the space Q£°, then B k is a left inverse provided its range is restricted to the same 
space. We introduce the notation A^) to denote a matrix A with its first k rows set to 
zero. Thus, we have that B k k ^D k = Iq oo . We note that these relationships carry over to 
finite truncations if we replace the last column of B and the last k columns of B k k ^ with 
zeros, since D k : Q k — ► Qo~ k while B k k j : Qo~ k — > It is easy to see that these 

simple inversion (integration) operators originate in the recursions 


( 3 ) 


rpt 
1 Jc + 1 

k + 1 


TL 1 

lb - 1 


= 2 T k , 


k = 1, - , 

, T{ = T 0 , 


rpn OT" T" 

(A\ ±k+2 ^±k , 1 k-2 rp r _ 1 

4(fe + 2)(fc + 1) 4(Jb 2 -l) + 4(Jb - l)(Jb - 2) “ * ’ 

r" = o , r" = o , r" = 4T 0 , 

and so on for higher derivatives. Clearly, B and B* 2 ) are the matrices of recursion 
coefficients for equations(3), (4) respectively. In the discussion we use the same symbol 
for an infinite dimensioned matrix operator and its finite dimensional truncation, where 
the distinction is clear from the context. 

More generally, if {£2fc(a0}£° is a family of orthogonal polynomials, then a three term 
recurrence for multiplication by the monomial x 

i 

(b) ) ] Qk+l&k+lik = %Qk j k = 0, 1, • • • 

l=-l 

follows easily from the orthogonality of the Q k [17]. Since the Q' k are orthogonal (with 
weight p(x), as is easily seen by integrating (2) by parts), they also satisfy a relation of 
form (5): 

(®) Q'k+l a i+l,k ~ X Q'k > * = 0,1,' •• 

I=-l 

Therefore, by differentiating (5) and combining with (6), we arrive at [7]: 


CO Qk+l^k+l,k — Qk ) k — 0, 1, 

/=-l 

which allows the efficient inversion of differentiation to all orders. The coefficients in (7) 
can be derived from those of the basic recurrence (5), which defines the f amil y 
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The method we shall present in Section 3, explained in detail for 2nd order operators 
but not limited to them, relies on restricting the domain of D n to the subspace Q n — 
spa 7 i{Qfc}^, thus ensuring the existence of a unique inverse. Throughout we tacitly 
assume that the operator L N , the N - th order Galerkin approximation to L, has rank 
N - n when acting on elements of Q$. Thus, the problem of solving the resulting 
algebraic system for right hand sides restricted to Qq n has a solution containing n free 
parameters. We moreover assume that the operator is nonsingular when restricted further 
to . Thus, the null space contains no element orthogonal to Qo These assumptions 
are not as restrictive as one might at first expect. The method is most effective when 
the above problem needs to be solved repeatedly for several right-hand sides / and high 
accuracy is desired. This type of problem arises, e.g. when the Navier-Stokes equations 
are solved in a geometry in which the Laplace operator is separable, and the boundary 
conditions are periodic in all directions except one. Common examples are provided by 
the Laplace operator in various separable curvilinear coordinates, where expansions of 
smooth functions in terms of eigenfunctions of the Laplacian in the bounded direction 
do not possess good convergence properties. 

In Sec. 3 we give some examples of the inversion of the Laplacian in some common 
geometries, including a disk and an annulus in cylindrical and helical coordinate systems. 
The use of the method for initial value problems is demonstrated through a study of the 
Airy equation, while the biharmonic equation analyzed in Sec. 4 provides an example for 
a higher order problem. Also considered is the Stokes problem: here a coupled system 
of second order equations is considered with boundary conditions given for only one. 
The method is easily extended to cover this case. The Chebyshev polynomials are an 
especially important family, because of their optimal approximation properties as well as 
the applicability of the Fast Fourier Transform. Thus, most of our explicit calculations 
are carried out for Chebyshev- Galerkin matrices. In Sec. 4 we carry out a detailed 
conditioning analysis for typical problems. It is found that if the leading coefficient m n (x) 
does not vanish in the interval under consideration, the method generically produces well 
conditioned operators. Finally, in Sec. 5 we discuss how to use rational mappings to 
stretch the coordinate system near points where the solution of a BVP exhibits rapid 
variation, thus ensuring a more efficient representation of the solution without sacrificing 
the speed of the method. 


2. Recursive determination of derivatives 

Throughout we assume that {Qfc(x)}^° is a family of orthogonal polynomials in [a, b\ 
with weight such that if u £ C°°[q,i b ] and if we set 

N 

UN = X) tfsQk 
0 
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with 

1 f b 

u° k = — u(x)Qk{x)w(x)dx , where hk =|| Qk ||^ 
fljc J a 

then the error || u — un lltu^-^ 00 0 faster than any negative power of N. This is for 
example true for the eigenfunctions of certain singular Sturm-Liouville problems [4]. 

We shall write for the restriction of the n-th derivative operator with respect to x 
on Qq = span {Qk}o ■ We adopt the notation 

(8) £>”<. = ■£>“ = 

0 0 

and we write ujv = col[u { ) € R N+1 , (t = 0, 1, ... ,N). In the sequel we will drop the 
subscript N when the distinction between truncated and untruncated expansions is clear. 
Also, as stated earlier, we will write A(k) f° r a matrix A whose first k rows have been set 
equal to zero. 

We now prove the following theorem, which is a special case of Theorem 2.2, but 
because of its simplicity serves to explain ideas. In this form, the theorem applies, e.g., 
to the Legendre polynomials: 


Theorem 2.1. If the family satisfies the recurrence 

(9) Qk+ i — Qk-i = f(tyQk > k = 0, 1, 

with Q-i = 0, then 


( 10 ) 

Proof: Clearly, 
so that 


M A 1 

*4-i 


f(k + 1) f(k - 1) 


= -tt® , k = 1,2, 


Q'i+i( x ) = Y f{m)Q m {x) 


u = 


m=0 

m+k even 


Y v-lQ'k = Y ulQk 

k= 0 k=0 

oo k—1 

= Ytfc 2 f ( m ) Q » 

k= 0 m=0 

m+J: odd 

oo I oo 

= Y \ f ( m ) Y *4 

m= 0 I fc=m+l 

V m+fc odd 


)Q» 
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and finally 

OO 

«m = f( m ) J2 

k—m+l 
k+m odd 

resulting in the recurrence claimed above. 

Applying the formula of Theorem 2.1 repeatedly, we can derive similar recursions for 
the inversion of higher derivatives as well. For example, for D 2 we have 


(n) 

.,(/(* +!) + /(* -!)) , 

/(* + l)/(Jfc + 2) Ui f(k+ 1 )/(*)/(* - 1) + f(k - l)f(k - 2) 


k = 2, 3, • • • 


The above formulae lead to simple algorithms for the computation of derivatives of 
functions expanded in terms of the Q’s as well as for the solution of simple I. or B.V.P.s. 
For example, the solution to the problem 


u x = 9{x) , ti(a) = a 


can be found, if 


in the form 


while 


<?( x ) = £ 5mQm(s), 


771—0 


Uu = 


9k - 1 


9k+l 


/(* - 1) /(* + 1) 


, i = l,2. ••• 


“0 = f“ - £ «•»$»(«)) /<?o(o)- 

\ m=l / 


Other simple linear B.V.P. of the form 

Lu = g , Bu = l 


ran be solved efficiently by the inversion of banded matrices if the D.O. L has constant 
coefficients. For example, let 


( 12 ) 


L = 


cP 

dx 2 


+ A 2 . 


In order to solve the B.V.P. (12) with B.C. 


(13) u(a) = a , u(6) = /? 

numerically, by assuming a truncated expansion for u(i) of order M, we set 

a 2 \ 2 *0 

\ Ufc — gk 
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in eq. (11) for k = 2, • • • , M to get, together with the r-conditions 

M 

f^u° m Q m (a) = a 

m= 0 
M 

m — 0 

an almost pentadiagonal system (except for the first two rows which are full) for the 
coefficients u This can be easily solved by LU decomposition. Thus the r-conditions 
are viewed as the first two equations, followed by the first M — 1 recurrence relations for 
the determination of the Uk , k = 2, • • • , M, with = 0. This is equivalent to the 

usual way of stating the T-method [10]. 

An alternative approach is suggested here. We specifically look for null vectors in the 
form ejt = Tfc+Ufc , Uk £ (Qo _1 )' L > & = 1> • • • > n— 1- Then, if u p € (Qo _1 ) x i s a particular 
solution, the solution to the BVP can be writen as u = u p + Y, a k£k, with a* satisfying a 
n x n system. Note that if repeated solution of the system is required with different RHS, 
the Uk need only be determined once, and there is a slight reduction in computational 
overhead of our method when compared, e.g., with an efficient implementation of the 
r-method. In fact, our method can effectively optimize the conditioning of a problem 
by restricting attention to the most stable subspace. So, for example, if one is required 
to solve the Poisson equation Au = —g in a region fi, where f2 is a 2- (3-)- dimensional 
rectangle with one (two) periodic directions and one bounded direction several times by 
adopting a Fourier- ( Fourier )- Q expansion the problem decomposes to equations of type 
(12) in the bounded direction for each Fourier mode. The LU decomposition can be 
performed in a preprocessing stage and the results stored, resulting in only ss (10M) 
operations per solution per Fourier mode at all subsequent stages. The cost is thus 
comparable to solving the Poisson equation in the pure Fourier case! Similar results can 
be easily derived for other O.D.O.s with constant coefficients. 

A straightforward generalization of the previous formulas, which is useful in deriving 
properties for the Chebyshev polynomials follows [7]: 

Theorem 2.2. If the family {Qjt(®)}^° satisfies the recurrence 

(14) jlQUi h w = Qk , * = o,i,-.. 

/=-i 

with Q-i = 0, then if f(x) = J2T-0 fkQk( x ) is a. sufficiently differentiable function 

(is) ? = Ew®,i=u- 

/=-l 

where the l-th derivative of the function f(x) has expansion coefficients fjfK 
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Thus, even in this more general case, the expansion coefficients of a function can be 
calculated from those of its derivative in 0(N) operations. The more general form in 
Theorem 2.2 will be useful dealing with Chebyshev polynomials, for which it agrees with 
the usual normalizations. The proof is straightforward, but we give it for completeness: 
Proof: We can introduce the vectors 

and 

g* = (<?0, <?!,•••) , <£ = (Qo, <?'!,•••)• 

Then, /(x) = q x fk and f = q' x fk- Also, by assumption, q' x B — q x where B is the 
coefficient matrix for recurrence (14). Combining we find 

(/ - Bp') = o. 

Assuming that the Q\ k \i = 1,2,- •• are independent (true for all families that satisfy 
eq.(5) ), we find the relation claimed. 

We note that the Chebyshev polynomials in the standard normalization satisfy (14) 
with bk,k ± l = (±l)/(2(Jfc ± 1)). Also the Jacobi polynomials in their standard normal- 
ization satisfy a relation of type (14). Strictly speaking, Theorem 2.1 applies only to 
the Legendre polynomials. (Although we can scale the symmetric Jacobi polynomials so 
that (9) applies). In any case, (14) shows that integration is always banded, and of a 
simple form (the recurrence coefficient matrix) for all the classical orthogonal polynomial 
families. The discussion following Theorem 2.1 was given to clarify ideas, and in principle 
could have been omitted. 

We now focus on the operator D. This operator has a one-dimensional null space, 
and if appropriately restricted, it has an inverse. An especially useful restriction involves 
the subspace Q f*. In this space, the operator D has a well defined inverse, which we 
will denote as B. Although D has a full upper-triangular matrix representation, B is 
banded. Indeed, assuming the recursion in the form of Theorem 2.2, we have that B is 
the coefficient matrix for the recurrence (14) (note that this matrix had zeroes in the 
first row since Q' 0 — 0). 

Similarly, D n must be restricted to Indeed, Jf(D n ) = Qq -1 , so that the operator 
D n is nonsingular on Q^, the orthogonal complement of its null space, and it has a unique 
inverse, denoted Bfo : Q%° — > Any two images of an element z £ under rc-fold 
integration differ by an element of Qo -1 . The specific form of fixes that element of 
Qq~ x to be the zero element. Thus, the solutions of 

Lu - ^(m k {x)D k )u - fix) , u € Q", 


( 16 ) 
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and 


( 17 ) 


Lu = £(m t (x)D l Bfa) z = /(*), 0 € C?“ 
fc = 0 


axe equivalent. Cleaxly, D k B^ ^ These operators differ since the second matrix 

has zeroes in its first n — k rows while the first has some non-zero elements there. 
Example: 


B h) 


and using 


operator B^ 

: Qo* -* 

Qf (for 

families that satisfy Theorem 2 . 1 ) is 

0 

0 

0 

. . . 

0 


0 

0 ^ 


0 

0 

0 

. . . 

0 

. . 4 

0 

0 


1 

Joh 

0 - 

L±h 

/s/2/1 

0 

1 

/s /4 


0 

0 


0 

1 

hh 

0 

h+A 

hhh 

0 

1 

hh 


0 

j 


0 

0 


0 

0 


0 

l 


fk-ifk-X 


0 

/fc-l~t~/*-t-l 
fk—X fkfk+X 


fk+xfk+2 / 


D = 


0 

f 0 

0 

fo 

0 

fo 


0 

0 

h 

0 

... fx 

0 


0 

. . . 

. . . 

. * . 

0 / 2 Jt-l 

0 


0 




0 

fak 

•••/ 


we find that 


#(i) / DB(2) — 


0 

/ 0 
0 

0 

0 


Ja- 

hh 

0 

x 

h 

0 


-fo 0 


- 4 - 0 


0 - 4 - 


0 

0 

0 


/fc -1 


/k+1 


) 


In the general case, the operator D is hard to write explicitly, as it is the ‘inverse’ (in 
the sense discussed above) of a general tridiagonal matrix. However, all that is needed in 
our method is the expression for D k B^ t which is identical to the matrix B*~ k k ) except 
for the first n — k rows which, in general, will contain some nonzero elements. These 
elements axe easy to compute however, as they can be expressed in terms of elements of 
B(j) and the n x n principal submatrix of D. For example, the operator DB ^ for the 
general case is identical to B(j) except for the first row, whose elements axe 

row o {dB{2 ) ) = doi (&io^ii> + &i2&2i> ^12(^11 + ^22)1 ^>12^23, 0 , . . . , 0 ) . 

Here, the elements for the classical orthogonal polynomials can be found in Table 
3 , together with the elements of the matrix A and other relevant quantities using the 
standard notation [ 1 ]. Also, doi is the corresponding entry of the differentiation matrix 
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D , which is simply the derivative of Qi expressed in terms of Qo- For example, for the 
general Jacobi polynomials, doi = a + /3 + 2, etc. 

The relations for the Gegenbauer polynomials can be constructed from those of 
the Jacobis since 

r (u) r (q + l)r(2a + n + l) p( a ,/?) 

r(a + n + l)r(2a + l) n 

where a = f 3 = u - 1/2. Arrays are indexed from 0 to N, the maximum order of 
truncation. 


3. The method 

We present now a method for the efficient inversion of operators resulting from the 
spectral solution of the ODE 


(18) 


Lu — ^2(mk(x)D k )u — f(x) , x £ = (a, &), 

k = o 


subject to the constraints 

Tu — c. 

The constraints are represented by a linear functional T of rank n. 

We assume now that the matrices M k , representing multiplication by Tn k , are banded. 
This is for example the case if the original D.E. had rational coefficients. After multiplying 
out the denominators, we are left with low order polynomial coefficients, and as a result 
of the simple recursion operator A of multiplication by the monomial x (5), these have 
banded representations as convolution operators. The simplest form is found if we expand 
the resulting polynomial coefficients in terms of the Qk, then exploit the properties of 
the banded operators of multiplication by Qk- lu constructing an approximate solution, 
we look for a solution of 

Lnun = fpj ; un £ Qo , /jv £ Qo ”• 

with Ln the Galerkin approximation to L, as usual [10]. The main result can be expressed 
as follows: 

Theorem 3.1. Assume that the M k are banded. Also assume that Q$ = Af(L N )(BQn ■ 
Then, if there is a solution un, it can be written as a combination of an element w £ 
Af(L N ) and an element u p £ Q ^ such that Lj^u p = /. The solution of the latter problem 
can be performed in O(N) operations . 

To construct the particular solution, let z — D n u p £ Q$ 71 so that u p = B( n ) z € Qn 
is uniquely defined. Then the equation can be rewritten 

= t, MkD'B^Z = F , 


( 19 ) 
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A A A 

where we have introduced U p , Z, F to represent the vectors of expansion coefficients. (We 
use the same notation, however, for the differential and integral operators as for their 
matrix representations.) Since a solution u p of this problem was guaranteed to exist, z, 
its n-th derivative, exists as well. Here we must note that in the case of weak solutions the 
highest derivative must be handled carefully, but in this case convergence would be slow 
and the method would be impractical. Now we address our main question: Is the new 
system any easier to solve than the original? This is clearly true: the integration operators 
are all banded, and to find u from z we perform one more banded matrix multiplication. 
To simplify the notation, in the rest of the paper, unless otherwise stated, we will write 
D 3 ~ n = D 3 j = 0 ... n — 1, where n is the order of the differential operator under 
investigation. 

Finally, we need to determine a convenient basis for the nullspace of the operator L\. 
We define 

e k = Qk + Wk , w k <E Q„, 

with 


Tivejt = 0 =>• LfjWk = —LnQic- 


Then 

n— 1 

Tu = Tu p + ^2 a k T ejt — c 
fc= 0 

so that when the numbers T k i — (T e k )i are found, we have 


r otkT k i = ci — ( Tu p )i . 
k = 0 

In other words, for every new RHS we simply need to solve the standard BVP for u p , 
then evaluate the quantities ( Tu p )i and solve an n x n system for the a k . The condition 
of this system is known in advance. 

Example: The radial Laplace equation for the n-th Fourier mode: 

+ (( s + °^) ~ n2u = f ; ue 

This leads to the pentadiagonal matrix 

(x + a) 2 1 -1- (x + a)D~ l — n 2 D~ 2 . 

Example: The helical Laplace equation for the n-th Fourier mode is 

This is transformed to the (almost nine-diagonal) matrix 

r 2 (l + r 2 a 2 )I + r(l - a 2 r 2 )D~ l — n 2 (l + a 2 r 2 )D~ 2 . 
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Example: The initial value problem for the Airy equation 

y - a 3 (x - x 0 )y = 0 , y{0) = .355029403792807 , y'(0) = .258819403792807a 

has the solution 

y(x) = Ai(a(x — Xo))- 

the Airy function of the first kind. Here we include the parameters a, x 0 G [-1, 1] in 
order to scale the interval over which the problem is solved, since Chebyshev expansions 
apply naturally over the interval [ — 1 , 1] . For x > 0 the solutions decay exponentially 
and the numerical algorithm converges rapidly. However, for x < 0 the solutions exhibit 
oscillatory behavior with ever increasing frequency, and convergence can only be achieved 
if sufficient modes are include to resolve the most rapid oscillations present. In Figure 1 
we show the solution to the problem with xq = — 1, a = 10 with N = 30,40 respectively. 
The first case is underresolved, and the maximum absolute error is O(10 -2 ), while the 
second case is barely resolved, and the error is 0(1O -5 ). A slight increase in the order 
of truncation improves the solution dramatically. With N = 64, the error is already less 
than 10 -11 . 

Example: The two-dimensional Stokes problem is expressed by the system 

(20) A = -u, 

( 21 ) Hu-f 

where Ajv is the non-penodic part of the Laplacian for the m-th Fourier mode in a two- 
dimensional geometry with one non-penodic and one penodic direction, likewise, H is a 
second order linear operator with rational coefficients. No conditions are given on u while 
^ and tp x are specified at x = ±1. Such a system results, e.g., from the time discretization 
of the Stokes equations in appropriate two-dimensional domains. We consider projections 
f f N _ A € Q%~ A , L J ->■ ujn-2 € Qo~ 2 and rf> -» G Qo • We determine a particular 
solution for (21) as u> p G Q 2 2 and homogeneous solutions uik = Qk "F flfcj k 0,1, with 
Q k G Q*-\ The general solution for (21) is then 

( 22 ) & N -2 = + <* 0^0 + 0 * 1 ^ 1 . 

The general solution of (20) can now be written as if ) n — tp p + fioTpo + AV’i where 
ij> k = Q k + <f> k+2 , k = 0, 1 (with ^k +2 G Q 2 ) are the homogeneous solutions and Vv = 
<5f p q. a 0 $ 0 + G Q 2 is a particular solution with A^ p = —<*> p and A = 

—cjk,k = 0,1. The boundary conditions can now be applied to to produce a 4 x 4 
system in the at, /3k, k = 0, 1: 

Act = </>, 

A u+ i = ^(l) , A 2 , j+1 = ^(-1) 

A 3J+ i = «*,(!) , A4j+1 = *j>(-1) 


with 
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(j = 1,2, 3, 4) so that A need only be evaluated once, ct* + i = a* , k = 0, 1 and &k + 3 = 
ftk , k = 0, 1 and <f>i = ^>(1) — 'fp(l) etc. If this problem is solved as part of a time 
integration of the Navier Stokes equations in a two-dimensional region, the only functions 
that need to be computed repeatedly are and u p G Q%~ 2 - 



Figure 1: Computed solutions with 30 modes (triangles) and 40 modes (circles) 
are plotted versus the exact solution of the Airy equation for a = 10. 

Finally we comment on the solution of the BVP with arbitrary BC. In this case, we 
need to determine the null space, and add an arbitrary combination of nullvectors to 
the particular solution to satisfy any desired BC. An alternative form of our method can 
also be considered. It is based on commuting the polynomials with the differential oper- 
ators and multiplying on the left by the integration operator D n , so that the differential 
operator matrix becomes banded. This approach is discussed in [7]. This is in fact the 
T-method, and various instances that have been worked out (eg. [10], [8]) are of this 
type. Theorem 2.2 establishes the success of this approach as a consequence of the basic 
recurrence relation (14). 

However, other preconditioners may be available, depending on the special structure 
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of the matrix operator L N [16]. The present method’s main appeal, besides the fact that 
a convergence analysis is available (see Sec. 4) is its simplicity and generality. Indeed, 
the complicated expressions for the differential operators are entirely avoided. Of course 
it should be expected that in special cases simpler forms and preconditioners might be 
possible. An example of this is provided by the Laplace equation in a circle for which the 
integration preconditioner leads to pentadiagonal forms while a simpler tridiagonal form 
is in fact possible [16]. This is related to the special properties of the operator (x + 1)^. 
We have not yet explored the flexibility of the choice of spaces, which may sometimes 
lead to more efficient algorithms. 

We must mention that the basic idea of the method presented here can be found m 
the monograph by Fox and Parker [9]. They discuss the solution of ordinary differential 
equations with low order polynomial coefficients by means of Chebyshev polynomial 
expansions, and they show that integrating the given equation as many times as its order 
leads to banded problems easily treatable by difference equation techniques. A variant of 
this idea, employing the backwards recurrence relations for the coefficients of a function 
in terms of the coefficients of its derivative, again for the Chebyshev polynomials, was 
used by Clenshaw [5]. 


4. Stability and convergence 

As differential operators are unbounded (in the usual norms), so are their difference 
and spectral analogues under mesh refinement. Numerical studies have shown that the 
spectral radii of Chebyshev and Legendre differentiation matrices are 0(N) where N is 
the dimension of the subspace. Moreover, these matrices are far from normal, so that 
their norms can grow even faster. For example, the maximum norm of the differenti- 
ation matrix, D, for a family satisfying the simple recursion (9) such as the Legendre 
polynomials, satisfies the lower bound: 

„ .N — m 

(23) \\D [|oo > max/(m) — - — . 

From Table 3 we see that f(m) = O(m) so that the lower bound above is 0(N 2 ). This 
lower bound grows accordingly with the order of the derivative being approximated. (For 
example the second derivative can behave as 0(N 4 ).) 

The poor conditioning of the matrices arising from spectral discretizations both limits 
the accuracy of solutions, due to roundoff errors, and imposes severe limitations on the 
time step for explicit solutions of dynamic problems [15] or on the convergence rate of 
iterative solvers. It is well known that the reformulation of differential equations as inte- 
gral equations often leads to bounded operators and well-conditioned problems. As our 
formulation of the discrete equations is based on integral operators, we also expect to ob- 
tain well-conditioned linear systems. For some constant coefficient problems, Greengard 
[11] has directly analyzed spectral approximations to equivalent integral equations and 
demonstrated the gains in accuracy which can be attained. In this section we generalize 
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and expand on Greengard’s results to cover the algorithms we’ve proposed. We use the 
stability estimates we obtain to prove convergence of the method. 

4.1. Estimates of the condition number. We assume the differential equation takes 
the form (18) where the coefficients rrij(x) are polynomials and mo is bounded away from 
zero on [a, 6]. Of course this last condition is required to avoid singularities in the solution, 
where the spectral approximation itself may not be well-behaved. We concentrate on the 
system (19), used to determine a particular solution: 


Mo + £ Z = AZ = F , 

where the Mj are the Galerkin approximations to multiplication by the polynomial co- 
efficients and D~ 3 — D n ~ 3 B* n y The additional problem to be solved, involving the 
boundary conditions, will be of much lower dimension, and its conditioning will depend 
on the specific constraint conditions. We begin by estimating various norms of the in- 
tegration operators. We view these, now, as operators on I 2 and loo- The bounds we 
derive obviously extend to finite truncations. We make the following assumptions about 
the orthogonal family, which are satisfied for proper normalizations of any of the Jacobi 
polynomials, as well as the Hermite polynomials (see Table 3): 



Assumption 4.1. The orthogonal family Qk satisfies: 


(25) 


su Pfc ^ Qk,Qk 2 ^ 

— — = Kq < OO 

inf* < Qk, Qk >w 


(26) P>°- 

Here, < , denotes the weighted inner product defining the family. The best exponent, 
p, is 1 for the Jacobi family and 1/2 for the Hermite family. 

(To use Table 3 to verify the statement concerning p, form the normalized families 

by dividing Qk by y/hk and note that bkj is transformed to bkj\J~hkjhj .) We have the 
following lemma describing the structure of the integration matrices: 

Lemma 4.1. The matrices D~* are banded with bandwidth j, with the possible exception 
of a finite number of elements in the first j rows, and there exists a constant Bj such 
that |(D“- , ’) W | < Bjk~ jp . 

Proof: We first note that the integration operator B^ = D~ n coincides with (B^)” 
except for the first n rows which are zero. Since B^j is tridiagonal with elements satisfying 

(26) , the result is immediate. For the other terms we use the fact that DB ^ X j = I to 
write: 

(27) D~ 3 = D n ~>Bl n) = D^HB^T +C) = (B^) 3 + D n ~ j C, 
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where the nonzero elements of C are simply the negatives of the nonzero elements of 
(B^) n in the first n rows. Since D n ~ j is upper triangular with nonzero elements only in 
superdiagonals n — j — ► oo we see that the nonzero elements of D^ n ~^C are restricted 
to the first j rows and 2n + 1 columns. Since the result holds for the lemma is 

proved. 

This leads immediately to the following theorem: 


Theorem 4.1. The operators D J : £ 2,00 h,<x> o- re compact. 


Proof: The boundedness of the infinity norm follows directly from Lemma 4.1 and the 
fact that the norm is equal to the maximum absolute row sum. For the 2-norm we have: 


\\D~‘yh = ( e £. 

(28) 


1/2 / 00 max(2n+l,t+j) 

<y/2^+lB } £ £ 

\t=l j=max(l,t- j) 

< (2n + l)Bj||y|| 2 . 


1/2 


W 2 


To prove compactness it is sufficient to show that D~ } can be approximated by a sequence 
of bounded operators of finite rank, { Df / }• For these we simply take the operators 
defined by setting all rows below M to zero. Repeating the arguments above we have: 


(29) || D~> - DJI H2.00 < (2 j + l)BjM~ jp -> 0 , M —> 00 , 


completing the proof. 

We would now like to bound the norms and condition numbers of the Galerkin poly- 
nomial multiplication matrices, M k - We begin by showing that the Galerkin matrix is 
nonsingular if the zeroes of m k lie outside [o, &]: 


Theorem 4.2. Let M k be the matrix representation of the Galerkin approximation to 
multiplication by the degree q polynomial m k {x) relative to the orthogonal system {Qj{x)}q 
on [a, 6]. If the zeroes of mk lie outside [a, 6], then Mk is nonsingular. 

Proof: Suppose the contrary. Then there exists a nonzero polynomial, P, of degree N such 
that mkP = YJj=n+i CkQk(x). We then have that m*.P is orthogonal to all polynomials of 
degree less than or equal to N and has at least q zeroes (counting multiplicities) outside 
[a, 6]. This implies that m k P has at most N zeroes of odd multiplicity in (a, b). Let r t , 
i = 1, . . . , s denote these zeroes. Then S(x) = n*=i( x - r «) is a polynomial of degree 
less than or equal to N such that m k PS is of one sign on [a, b). However, we also have 
/ a 6 Ljm k PSdx = 0 by the orthogonality of m k P to polynomials of degree not more than 
N. This is a contradiction, so P cannot exist. 

An immediate corollary of this theorem is: 


Corollary 4.1. The spectrum of M k is contained within {y = m^i), x € [a, &]}• 
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Proof: Suppose Mk — XI is singular. Since M* — XI is the Galerkin approximation to 
multiplication by m* — A, we conclude m* — A must have a zero in [a, 6], completing the 
proof. 

From the eigenvalues we can easily bound the norms: 

Theorem 4.3. The matrices, Mk, satisfy the following bounds: 

(30) \\Mkh < «<? max |m*(x)|, ||M fc -1 || 2 < k q max |(l/m*(x))|. 

xe[a,b] xG[o,6J 

Proof: Let Mk denote the matrix representing the Galerkin approximation to multipli- 
cation by mk relative to the orthonormal basis obtained by normalizing the Q' s. Since 
Mk is symmetric, its 2-norm and the 2-norm of its inverse are bounded, respectively, by 
the largest and the inverse of the smallest eigenvalues (in absolute value). These are 
in turn bounded by max x6 [ ai i,] |mfc(x)| and maXjg^j] |(l/m/t(x))| from Corollary 4.1. Let 
R = diag(^/< >u,)- Then Mk = R~ x MkR. Taking norms yields the final result. 

We have shown that the system defining the particular solution has the form Mo + K 
where, for regular problems, Mo has a bounded condition number uniformly in N and K 
approaches a compact operator. We have, then, the following alternative: 

Theorem 4.4. Either there exist constants Co and C\ and an integer No such that for 
all N ^ Nq and vectors, y, wxth T/Uchdean norm 1, 

(31) Co < \\Ay\\ 2 < C x 

or 0 is an eigenvalue of finite multiplicity for the limiting system. 

Proof. Suppose that zero is not an eigenvalue of the limi ting system. This system can be 
written in the form Mq(I + K) where K is compact and M 0 is bounded with a bounded 
inverse. Then zero is not an eigenvalue of I-\-K, which must then have a bounded inverse. 
For A we have A = Mo(I + Kn) where Kn, when viewed as an operator on / 2 , satisfies 
\\K — Kn\\ —* 0, N — > oo. By the Banach lemma (I—Kn) = (I - K +(K — Kj*)) will have 
a uniformly bounded inverse for N sufficiently large. It is also uniformly bounded above. 
Since we have uniform upper and lower bounds on Mo, the existence of Co,i immediately 
follows. 

We remark that the zero eigenvalue will occur only if there exists a nontrivial solution 
of the homogeneous differential equation which is orthogonal to all polynomials of degree 
less than n. This is clearly not generic, and if it holds it can be remedied by looking for 
particular solutions in a different subspace. In the future we plan to consider the case of 
singular limiting problems in more detail, particularly in cases where the lead coefficient 
is zero somewhere in [a, b ]. 
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4.2. Error estimates. Given these bounds on the condition number of the linear sys- 
tem, a convergence result is easily proved. We explicitly assume that the original problem 
has the following properties: 

Assumption 4.2. 

(a) : The constraint/boundary operators T satisfy an inequality of the form \Tw\ < 

IMU- 

(b) : The forcing function, f(x'), is in C T {\a, 6]), r > 1. 

(c) : If w is a solution of the homogeneous problem ((IS) with f — 0) satisfying 
Tw = 0 or <w,Qk >w= 0 for all k = 0, . . . ,n — 1 then w = 0. 

We now prove a sequence of estimates of various parts of the error. For the continuous 
problem, Assumption 4.2 implies the following (e.g. [6]): 

(i) : There exists a basis, {'u J }"=o € C^Qo, b)), for the space of solutions to the 
homogeneous problem taking the form Uj = Qj + Uj, with Uj orthogonal to Q£ -1 . 

(ii) : There exists a unique solution, u € C n + T {[a, b]), which can be written u(i) = 
u s (x) + 5Zj=o c j u j( x ) u s orthogonal to Qq~ 1 ■ 

(iii) : The n x n matrix 

T e — T uq T u \ ... T tt n _i j 

is nonsingular. 

Let u 4 (x) denote the approximate particular solution, that is, the polynomial whose 
expansion coefficients are given by D~ 2 z , and Vj(x) denote the approximate solution of 
the homogeneous problem taking the form Qj+Vj with Vj orthogonal to Qo~ l ■ We assume 
that the right-hand side of the inhomogeneous equation is obtained via interpolation at 
the relevant Gauss, Gauss- Radau or Gauss-Lobatto points. Let 

(32) T a = Tv o Tv i ... Tv n - X ] . 

We then have: 

Lemma 4.2. There exists No such that for N > Nq: 

(i) : There exist constants Gi t(i such that: 

\\u { P - t>f||u, < G^N-*, 0 < l, fi < oo, j = 0,... ,n- 1. 

(ii) : There exist constants R ^ such that: 

\\T e -T a \\ <R»N-\ 0 < /z < oo. 

(iii) : There exist constants D\ such that: 

Ik" - k’lk < AJV>-'||/|U,„ o < i < n. 
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Proof. We rely extensively on the approximation results listed in [4, Ch. 9]. Now 
Uj — Vj = Uj - Vj. Let Uj = Uj t N + Wj where ufy G Qo~ n and G Q^-n+i ■ Then: 

(33) uj - vj = (uj,N - Vj) + xbj. 

Estimates of the last term and its derivatives follow directly from results on approximation 
by singular Sturm- Liouville eigenfunctions and the smoothness of Uj. For the first, we 
rewrite the expansion coefficients as B^Z u j 7 n and B^Z v j and introduce Z ej = Z Ut j^ — 

Z V j . Let Z Ut j t jy = Z Ut j,N+Q — En+qZuJ,n where Q is the bandwidth of the matrices A 
and E represents extension by 0 of a vector to a longer vector. Denoting explicitly by 
Am the matrix A associated with degree m + n — 1 truncations and by P m the restriction 
of a vector of order larger than m to the m-vector containing its first m components we 
have: 

(34) As-nZe ,j = Pv- n AN- n+ QZ u j'K. 

By the properties of A we have: 

(35) \\tjh < C\\Z^h- 

Now Z u j t w can be estimated by derivatives of Wj. Therefore, we have estimates of the 
nth derivative of the error in terms of the difference between u) ’ and its projection into 
Qo~ n . From [4] we directly obtain the estimate in (i). For lower derivatives we simply 

A 

apply the bounded operators to For higher derivatives we apply the derivative 
operators, which, though unbounded, still contribute only polynomial growth. To derive 
(ii) we use the estimates in (i) and assumption (a) on the constraint operators. 

Estimates of the particular solution follow the same pattern. Introduce /y, the poly- 
nomial approximation to / used to compute v S} and /y, the orthogonal projection of / 
into Qo~ n . Write u[ n ^ = where G Qo~ n and € Q^-n+i- Let the 

expansion coefficients of u 3} n and v 3 be given by B^Z SjUi n and B^Z SiV respectively. Let 
R 3 = Z s v Z 3U and Z 3M ^ r = Z s u j^^.q E>j ^.q Z 3 ^ u j ^ t . Then we have: 

(36) AN- n R 3 = PN- n AN-n + QZ 3tU + Fjv — Ftf. 

From the boundedness of the .A’s, the first term can be estimated in terms of ||u4 n *||w = 
0(W -r )*||'u>£ n )|| av .. (This holds because tx S) jv was constructed by projecting into Qo~ n 
and applying integration operators.) The second is bounded by the sum of ||/ — /y|| w = 
0(N~ r ) • ||/||u,,r and \\f — JnWu = 0(Ni~ r ) • The boundedness of A~ l yields then 

implies (iii) for the nth derivative. The bounds for the lower derivatives then follow by 
application of the bounded integration operators. 

We are now in a position to prove: 
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a = 5 

a = 10 

a = 20 

N 

<?l 

0AT-1 

k 2 

C7 

&N - 1 

*2 

<?i 

O’N- 1 

«2 

32 

46.3 

.077 

605 

374 

.007 

53517 

2992 

.100 

29891 

64 

46.3 

.077 

605 

374 

.023 

16015 

2992 

.042 

71295 

128 

46.3 

.077 

605 

374 

.023 

16015 

2992 

.008 

378611 

256 

46.3 

.077 

605 

374 

.023 

16015 

2992 

.008 

378611 

512 

46.3 

.077 

605 

374 

.023 

16015” 

2992 

.008 

378611 

1024 

46.3 

.077 

-4 1 

605 

374 

.023 
1 1 

16015 
1 

2992 

■ . ~ 3f 

.0081 
_ 1 1 \ 7 

378611 

-^—2 


Theorem 4.5. For some No < oo there exist constants Hi such that, for all N > No, 
the difference between the true solution, u, and the approximate solution, v, satisfies. 

||u (0 - v (,) |U < HiN*- r \\f\U,r, l = o, . . . ,n. 

Proof. We have, 

n— 1 

(37) u = u, + Y, T e 7 = c - T u s , 

3=0 
n— 1 

(38) V = v s + T a 6 = c- Tv s , 

j=0 

Introducing e = u — v, v = / y — 6 and talcing the difference of the equations above we 
obtain: 

(39) e = u s — v 3 + ^2 (7 j{uj — vj ) — VjVj ) , T a v = ( T a — T e ) 7 — T(u s — v s ). 

i = 0 

Applying estimates (ii) and (iii) of Lemma 4.2 and the Banach lemma to the second 
equation we obtain \v\ = 0(iV2 -r )||/|| u ,, r . Substituting this into the first equation and 
again using parts (i) and (iii) of lemma 4.2, we obtain the desired result. 

We note that the estimate for the nth derivative is of optimal order for finite r. Of 
course, for / £ C°°([a, b]), we have convergence at a rate faster than any negative power 
of N. 

4.3. Direct computations of the condition number. Finally, we illustrate the con- 
ditioning results by computing the singular values of the matrix used in the numerical 
example in Section 3, namely Airy’s equation with a Chebyshev discretization: 

(40) A = I + a 3 { x + 1)D- 2 . 

The singular values for various N and a were computed using the lapack routine, dgesvd. 
The results are presented in the table below. 
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a = 1 

Q 

II 

k- 4 

O 

O 

a = 10000 

N 

<Tl 

&N - 1 

k 2 


<7JV - 1 

k 2 




32 

1.00 

.995 

1.01 

1.31 

.602 

2.17 

69.9 

.070 

1004 

64 

1.00 

■EEE1 


1.31 

.602 

2.17 

69.9 

.070 

1004 

128 

1.00 

.995 

1.01 

1.31 

.602 

2.17 

69.9 

.070 

1004 

256 

1.00 

.995 

1.01 

1.31 

.602 

2.17 

69.9 

.070 

1004 

512 

1.00 

.995 

1.01 

1.31 




.070 

1004 

1024 

1.00 

.995 

1.01 

1.31 




.070 

1004 


TABLE 2. Extreme singular values for I — aD 4 . 


We see that the extreme singular values and, hence, the condition number of the system 
matrix are independent of N, once N is taken large enough to resolve the problem. 
(The large condition number for large a simply reflects the large but bounded condition 
number of the integral equation.) To illustrate the insensitivity of this result to the order 
of the underlying differential equation, we have carried out the same computation for the 
biharmonic; that is for A = I — aD ~ 4 , with the results tabulated below. 

Here, the results are quite independent of the truncation, as the extreme singular values 
are resolved with N = 32, so the only growth in the condition number is associated with 
the growth of a. 


5. Rational Maps for Layer Resolution 


The Chebyshev approximation to a function with a region or regions of very rapid 
variation may exhibit Gibbs-type phenomena, that is large amplitude oscillations of the 
error, unless many basis functions are used. Therefore, adaptive computations using coor- 
dinate mappings to stretch these regions have been proposed. Bayliss and Turkel [3] have 
made a comparative study of various functional forms, all of which were transcendental 
functions. 

In order to be able to use our fast solvers, we consider rational maps. That is, we 
directly solve ( 1 ) in y space where: 

_ P(y; 

the polynomials P and Q are of low degree, and 77 is a parameter vector. In an adaptive 
procedure, 77 would be chosen to minimize some measure of the error, for example the 
error functional proposed by Bayliss and coworkers [2]. A very simple construction of 
an appropriate map can be motivated in the following way: Let g(y) = P/Q ■ The 
convergence of the Chebyshev expansion (in y ) depends on the behavior of 



d k u 

dy k 



( 42 ) 
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This suggests that improved convergence will follow from making dg/dy small where 
d k u/dx k is large. We imagine an underlying linear map (so that the limits of the com- 
putational region will be [ — 1, 1]) stretched near a finite number of points, Xj. This can 
be accomplished by subtracting scaled and shifted multiples of the function: 

ay + 7 

(43) h,(y\ a, (3, 7 ) = j 

Note that h' s (0] a,/3,7) = a and that the derivative approaches zero as (3y 2 — * 00. We 
then propose: 

(44) g(y) = Sy + C — ^ h»(v ~ VA °y> #>> Ti) 

i 

The n umb er of terms in the sum, and, hence, the degree of the map and bandwidth of 
the resulting matrices, depends on the number of layers present. Then, if S - ctj is small 
and j3j is large, enhanced resolution at g{yj) will be obtained. 

We demonstrate the idea on the following simple boundary value problem: 


(45) + X Tx = °’ ^ (±1) = ±1 - 

This problem has the exact solution 

#(x) = -1 + 

which, for e small, exhibits a region of rapid transition near x = 0 whose width is l 

We see below that for e small and even a large number of modes, there is a very strong 
Gibbs-like behavior. We consider the rational map 


(46) 


2 A + y 2 
X ~ A+1 V 1+y 2 ’ 


which is derived from the general expression above, making use of the symmetry. In 
particular, the derivative is minimized at y = 0 where its value is 2A/(A + 1). Under the 
change of variables, ip{y) = <£(i(y)), the equation becomes: 


(47) 


(Pip 

6 W + 


«*>■ -<$<>) Ty =“■ 


± 1 . 


We studied this equation for various parameter values. We found that with e = 10 -12 , 
a value of the parameter A of the map of the order of 10 -6 yielded the best results. This is 
reasonable, as one might expect A = 0(y/e) to match the scaling in the layer. We did not 
systematically search for the optimal value. In Figure 2 we present the solutions obtained 
for various numbers of modes, N , and mapping parameters, A. Note that A = 1 yields 
the identity map, i.e. the case of standard Chebyshev approximation. With A — 1 and 
N = 256 we see oscillations near the layer with an overshoot of about 18%. Increasing 
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N to 32768, which was the largest value considered, had no effect on the amplitude of 
the overshoot, but did contract the region of oscillation. With A ^ 1, on the other hand, 
we obtain reasonable results with many fewer nodes. For example, with A = 10~ 6 and 
N = 64 the overshoot is less than 1%. Increasing N to 1024 and spending a bit more 
effort optimizing the parameter reduces this to 3 x 10 -4 . 


I I 



I I I l _ _l I I I l _ 


Figure 2: Solutions for e = 10 12 . 

It is clear, that there is no change in the essential (O(N)) amount of work needed to 
solve the problem, although the bandwidth does increase by approximately a factor of 
5. If an iteration were used for the minimization of some error functional, by shifting 
the position of the shock and changing the magnification factors, each step would require 
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the recomputation of the operator coefficients, and the solution of the problem. These 
procedures are of comparable numerical cost, so that the desirable features of the method 
axe essentially preserved under the change of variables. 

It is worth noting that there are limits to the capability of this method to concentrate 
a large fraction of the mesh in a small region. A simple calculation indicates that g (y) = 
0(e) in a y-interval of width y/e. To achieve a greater magnification one must use rational 
maps of higher degree, which results in larger bandwidths. A more detailed study of the 
properties of rational coordinate mappings is planned for the future. 
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